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I. INTRODUCTION 

The properties of a variety of astrophysical objects are 
essentially determined by the high-density equation of 
state (EOS) of hydrogen^. Hydrogen is fully ionized 
at these high densities independent of the temperature. 
In this metallic phase, the electrons can be of arbitrary 
degeneracy including the highly degenerate limit (T = 0) 
and the ions are strongly coupled. The thermodynamics 
is thus strongly influenced by the well-pronounced short- 
range order in the proton subsystem although the major 
contributions stem from the electron gas. 

Metallic hydrogen occurs in the interior of giant 
gas planets like Jupiter, Saturn and similar cxtrasolar 
planetai^. The temperatures encountered along the 
isentrope of giant planets are in the order of a few elec- 
tron volts. At densities comparable to solids and above, 
a T = description of the electrons is possible for colder 
planets, whereas temperature related corrections might 
be needed for hotter planets^. Higher temperatures and 
densities are found in white dwarf stars and the crust of 
neutron stars^. As most elements are fully ionized un- 
der these conditions, matter behaves hydrogen-like and 
the EOS is again determined by a combination of degen- 
erate electrons and strongly coupled ions. 

Beyond astrophysical applications, the EOS of dense 
hydrogen is also required to model inertial confinement 
fusion experiments as the compression path of the fuel 
runs through the parameter space considered here^. Al- 
though at much higher densities, the fully compressed 
DT-pellet has similar properties when considering elec- 
tron degeneracy and ionic coupling strength. 

Most astrophysical objects exhibit an isentrope that 
covers several phases. For example the internal structure 
of giant gas planets is determined by phase transitions as 
the molecular to metallic or atomic to metallic transition 
in dense fluid hydrogen (sometimes named plasma phase 



transitio n 10 ' 11 ). After many investigations leaving the 
nature of this transition unclear— ~—, recent first prin- 
ciple simulations showed that a first order phase tran- 
sition with surprisingly small volume change is indeed 
likely^. 

However, present ab initio simulations like DFT-MD, 
path integral Monte Carlo (PIMC), or coupled electron- 
ion Monte Carlo (CEIMC) cover a limited parameter 
space only. For a consistent description of the EOS, one 
has to require that i) different techniques agree in over- 
lapping regions and ii) the simulation data merge with 
well-founded theories in limiting cases. The second point 
has been achieved only in the low density region where 
PIMC simulations match density and fugacity expansions 
of the EOS perfectly^. 

In this paper, we focus on the high-density limit of the 
equation of state of fluid hydrogen. As quantum Monte 
Carlo schemes are not available for very high densities, we 
rely on DFT-MD simulations here. However, first prin- 
ciple simulations like DFT-MD have so far been unable 
to provide results that converge into the high density, 
i.e. Gcll-Mann & Brucckncr, T — limit. We resolve 
this issue by carefully performed DFT-MD simulations 
for higher densities and by employing an analytic ap- 
proach that extends the T = limiting law to parame- 
ters with finite temperatures. We can then demonstrate 
agreement in overlapping regions of density, so that the 
goal of a combined EOS of hydrogen, that is solely based 
on methods in the physical picture and covers the entire 
density range for temperatures of a few electron volts, is 
reached. 

The analytic EOS theory we apply for high densities 
is a two-fluid model based on a perturbation with re- 
spect to the electron-ion interaction. It keeps all con- 
tributions from correlations in the ion subsystem and 
is applicable for arbitrary degeneracy of the electrons. 
Thus, our two-fluid model is valid for fully ionized plas- 
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mas with, compared to thermal energies and/or ion- ion 
interactions, weak electron-ion interactions. 

After an introduction of the model in the next section, 
we present results and comparisons with data from first 
principle quantum simulations. In particular, we show 
which steps are necessary to reach agreement between 
our two-fluid model and the simulations and also give 
limits for the applicability of both. 



II. ANALYTICAL EOS APPROACH 

We consider fully ionized plasmas consisting of protons 
and electrons. To characterize the interaction strength 
between the particles, we define the classical coupling 
parameter 



with 



r = 
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In the quantum regime, the classical kinetic energy scale 
has to be replaced by its quantum analog: fcgT — > | (K a ) . 
The mean kinetic energy (K a ) can be calculated via a 
Fermi integral 2 ^ which recovers both the classical as well 
as the fully degenerate (Fermi energy) limits. Note that 
classically all coupling parameters arc identical for hy- 
drogen while the electron-electron and the electron-ion 
coupling are strongly reduced in plasmas with highly de- 
generate electrons as (K e ) ^> fc^T holds here. In the 
quantum limit, the electron coupling parameter T e be- 
comes smaller than unity due to quantum degeneracy. 
Here, the Brucckner parameter r s = d/as is commonly 
used to describe the coupling strength. 



Two-Fluid Model 



F = F eg + -1 



Ni 
2 



Vei(k) Xee(k) 



(3) 



A similar model was also used by, e.g., Chabricr & 
Potckhi n 28 ! 29 . The first term in Eq. ([2]) represents the 
ideal contribution from the classical ions. In the second 
term, the density derivative of the free energy F pro- 
duces the contribution of the correlated electron gas via 
the free energy of the isolated electron subsystem F eg . 
The integral over the electron density response function 
Xee, to be taken in random phase approximation (RPA), 
is an electron-ion cross term that arises from the linear 
response treatment of the electrons in the two-fluid de- 
scription. The third term in Eq. ^ accounts for ionic 
correlations described by the pair distribution gu. The 
polarizability of the electron gas is here to be taken into 
account via an effective ion- ion potential vff . This po- 
tential must also be applied when calculating the ionic 
pair distribution. The density derivative arises as the 
effective ion-ion potential is density dependent via the 
screening function. 

The internal energy can be calculated similarly via 
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Again, we have first the ideal ion contribution and C/o 
denotes the internal energy of the electron subsystem. It 
may be approximated by the free energy F ^ for highly 
degenerate systems near T = 0; otherwise it is given by 
Uq = F — TdF/dT. Ionic correlations are accounted 
for by the integral term. The additional temperature 
derivative is due to the fact that the effective ion-ion 
potential is also temperature-dependent. 



For the parameters considered here, the electron- 
proton interactions are weak (r e i<l) while the proton- 
proton coupling is usually strong (Tu > 1). The weak 
interactions between electrons and protons allows us to 
apply a Born-Oppcnheimer approximation and treat the 
electrons in linear response to the fields created by the 
ions. As a result, one can rewrite the full two-component 
problem as a one-component system for the ions which 
interact via effective potentials. More precisely, this will 
be called the two-fluid model where the fully correlated 
electron and ion fluids interact only weakly with each 
other. Of course, this procedure eliminates the ability 
to describe bound states which is however unnecessary 
in the high density limit. Applying this two-fluid model, 
the pressure is given b y 26 > 27 

Pp = i | p dF(m) 

rii V drii 

- i^/*to,(D-i](r|- B4 A)#(r). 



B. Properties of the Ion Subsystem 

The effective ion-ion interaction consistent with the 
model above is given by 

vf (k) = v tl (k) + [v el (k)] 2 Xee (k), 

= ^e\k). (5) 

Here, the electron part in the static dielectric function in 
RPA, e^ 1 = l + u e eXee, was introduced. v ee = Aire 2 /k 2 
is the Coulomb potential between electrons 2 ^. The bare 
Coulomb interaction between the protons Vu is thus lin- 
early screened by the electrons. More simple approxi- 
mations for the effective potential can be obtained for 
small k, where the Debye or Yukawa potential, vff(k) = 
47re 2 /(fc 2 + k 2 ) with k 2 = (4e 2 m e /nh 3 ) dpf E (p), fol- 
lows. 

The derivation of the two-fluid description ([2]) and (|4]) 
clearly shows that the ionic pair distribution gu must be 
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obtained from a one-component description of the ions 
where the forces are given by the effective interaction i>|P. 
Possible methods to determine gu are classical Monte 
Carlo and molecular dynamics simulations or techniques 
based on integral equations like the hypernetted chain 
equations (HNC)^— . 



C. Contributions of the Electron Gas 

To describe the electron gas, we employ the quan- 
tum statistical method of thermodynamic Green's 
function o 25 ' 36 . Its advantage is the ability to describe 
systems with arbitrary temperatures including the cor- 
rect T = physics, the transition to Boltzmann statis- 
tics, and the correct high temperature (Debye-Hiickel) 
law. Using this technique, a perturbation expansion in 
the interaction strength can be establishe d 23 ' 25 . Includ- 
ing terms up to the second order, one obtains 

Pee(T e ,fI e ) = P l e d (T e ^ e )+p? F (T e ,f, e ) 

+pH w {T e ^ e )+pt n {T e ^ e ). (6) 

The terms arc the ideal gas law, the Hartree-Fock 
(HF) quantum exchange term, the direct Montroll-Ward 
(MW) term, and quantum exchange contributions of the 
second order (e 4 n), respectively As this perturbation 
expansion exists in the grand canonical ensemble, the 
density n e is related to the chemical potential fi e via 

n e (T e ,n e ) = — — . (7) 

An inversion within golden rule is performed in order to 
obtain the pressure as function of density^ 3 -. This means 
that correlation contributions of the free energy as func- 
tion of density are taken to be equal to the negative excess 
pressure as function of the chemical potential. 
The ideal pressure is given by 



2k B T ( 



pl d (T e ,He) = A 3 E h/2(Ve/ k B T e ) . 



(8) 



Here, A e = \j2ith 2 / m e kBT e is the thermal deBroglic 
wavelength and J3/2 is the Fermi integral of order 3/3^. 

First order exchange contributions are contained in the 
HF term 
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which is given by an integral over a Fermi integral of the 
order -1/2. 

The Montroll-Ward term can be computed by a double 
integral over the dielectric function of the electron gas e e 

00 00 
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It is consistent with the expansion ((6j) to use the RPA 
dielectric function e e . 

The normal e 4 exchange term, accounting for exchange 
effects of second order, can be written as an integral over 
Fermi functions, f p = [exp(/3p 2 /2m e — f3fi e ) + and 
Pauli-blocking factors, defined by /p = [l— fp], 



p% n (T e , Me ) = m, 
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Here, v ee is the bare electron-electron Coulomb potential. 

The expansion ([6]) accounts for direct correlations and 
dynamic screening, incorporates collective oscillations 
(plasmons) as well as quantum diffraction and exchange 
in the electron subsystem. This expression is valid for 
weakly coupled electrons of arbitrary degeneracy and in- 
cludes in particular the low and high temperature limit- 
ing cases of Debye-Hiickel and Cell-Mann & Brueckncr, 
respectively^ 3 -. 

Within the same quantum approach, protons can be in- 
corporated and an EOS for hydrogen can be calculated^ 3 -. 
The advantage of an EOS of hydrogen fully based on 
Green's functions is its capability to describe quantum 
effects in the proton subsystem correctly. However, a hy- 
drogen EOS based on expansion ((SJ) is restricted to weak 
coupling in the proton subsystem as well, a limitation 
avoided within the two-fluid approach presented here. 



D. Limiting Behavior and Applicability 

The condition of weak electron-ion coupling as used to 
derive Eqs. ((2|) and (j4]) is fulfilled not only for the high 
density limit of highly degenerate electrons and strongly 
coupled ions, but also in the high temperature and low 
density limits. Here, all interactions are weak and of the 
same order and the EOS is given by the Debye-Hiickel 
lawS 5 - 



Pp^Pp^ + (3p^ = J2n a -^- 



(12) 



The first term is the ideal classical gas contribution, the 
second term is the Debye-Hiickel correction determined 
by the total inverse screening length k 2 = ^ a 47rZ a e a n a /3. 
The sums in Eq. (fT2")) and in the definition of k run over 
all species a = {e, i}. 

This limiting law is not fully reached by our two-fluid 
model. However, it deviates only by a tiny fraction of 
p DH / p 2- fluid = 16^/2/23 = 0.99 (similar for the internal 
energy and other thermodynamic functions). The slight 
disagreement can be traced back to the neglect of the 
influence of the ions on the electrons. Note that this 
result can only be obtained if the contributions of the 
electron gas are evaluated for finite temperatures and 
not just in the ground state. 
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FIG. 1: (Color online) Lower panel: pressure of hydrogen 
as predicted by the two-fluid model ([2]) normalized by the 
ideal contribution. The high temperature (Debye-Hiickel - 
DH) and high density (T — 0) limits are given for comparison. 
Upper panel: ratio of the two-fluid model and the quantum 
statistical perturbation theory using the Montroll-Ward ap- 
proximation for hydroge n 23 ' 25 . 



Our quantum treatment via thermodynamic Green's 
functions also ensures that the electron contribution 
reaches the correct high density T = limit. Moreover, 
the electron contribution is the by far largest term in 
the thermodynamic functions for high densities. Thus, 
the two-fluid model also recovers the high density T = 
limit (with an error given by the ratio of electron to ion 
mass). 

Consequently, the two-fluid model constitutes a valid 
approximation to the EOS of fully ionized and weakly 
coupled hydrogen plasmas with arbitrarily degenerate 
electrons. Both conditions are fulfilled for temperatures 
above T = 2.5 x 10 5 K in the entire density range, where 
bound states do not occur and the coupling is sufficiently 
weak for the entire density range. 



E. Results of the Two-Fluid Model 

Figure [1] gives a general overview of the high-density 
EOS of hydrogen as calculated within the two-fluid 
model. Due to the normalization by the ideal pressure, 
the lines all approach unity for very high densities where 
the pressure is also independent of the temperature. At 
intermediate densities, the correlations yield a reduction 



compared to the ideal pressure. These correlations can 
result from the occurrence of bound states or occur be- 
tween free particles. Bound states can be excluded for the 
two highest temperatures in Fig. [1] Still, the improved 
description of ion-ion correlations within the two-fluid 
model yields a 10% correction compared to a description 
of the hydrogen EOS based entirely on Green's function 
theory. 

The two-fluid model of Potckhin & Chabrier gives al- 
most identical results to our approach in the parameter 
range where such a description can be expected to be 
valid, i.e. for metallic hydrogen and the high tempera- 
ture low density case^l. This is a nontrivial finding as 
the electron theories used in this paper and by Potckhin 
& Chabrier are quite different. 

In the case of the Green's function EOS 23 ' 25 , the cor- 
rect high-density limit is an intrinsic feature related to 
the quantum treatment of both electrons and ions. In 
the two-fluid model, only the electrons are treated fully 
quantum statistically. Still, the two fluid system shows 
a similar behavior as the high density limit is dominated 
by the electron contribution. This fact holds although 
the ion-ion correlations strongly increases with density 
since the electron terms grow significantly faster. In ad- 
dition, this behavior also minimizes the importance of 
the partition function of the strongly correlated ion fluid. 
For instance, small differences between Monte Carlo and 
HNC treatment of the ions at Ta = 130 give a deviation 
of 2% in the ion system which will change the result of 
the total equation of state by 0.16%. 

The approach to the exact high temperature limit 
(here given by the Montroll-Ward approximation of the 
Green's function theory) is demonstrated in the upper 
panel of Fig. [T] The two-fluid model can reproduce the 
exact law within 1% for T = 2 x 10 6 K. The essential in- 
gredient to achieve this result is the finite temperature 
description of the electrons. Moreover, it shows that the 
neglected influence of the ions on the electrons is very 
weak. 



III. DFT-MD SIMULATIONS 

DFT-MD simulations are a well-suited technique for 
the description of metallic hydrogen as it is encountered 
for example in the inner regions of giant gas planets. So 
far, there has been no overlap between an EOS calculated 
by DFT-MD and the correct T = high-density models. 
To close this gap and to allow for a direct comparison 
of first principle DFT-MD simulations to our two-fluid 
model, we carried out DFT-MD calculations for very high 
densities which require special adjustments. 

We used the DFT-MD programs VASP, CPMD, and 
abinh—ir— . The number of electrons and protons in the 
simulations was N = 128 . . . 432. The temperature of the 
protons was controlled by a Nose-Hoover thermostatic. 
The time step in the MD simulation was chosen to be 
At = 8a.u. = 0.194 fs. Every run covers at least 2ps 
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FIG. 2: (Color Online) Proton-proton pair distributions in hy- 
drogen at a temperature of T = 10 4 K for three different den- 
sities. Results of VASP (black dashed), CPMD (blue dashed, 
long spaces) and HNC using a linearly screened Yukawa po- 
tential (red solid) are compared. The black vertical lines in- 
dicate half the box length for the VASP runs. The lines for 
r s — 0.5 and r a = 0.7 have been shifted along the y-axis to 
improve readability. 

after an initial equilibration. The exchange correlation 
functional was of Pcrdcv-Burke-Ernzerhof type 4 ^. 

In VASP runs, the electron-ion pseudopotential was 
the standard (hard) projector augmented wave (PAW) 
potential as provided with the packag o 46 ' 47 and used here 
with a plane wave cutoff of 35 Ha. For densities higher 
than n = 2.2x 10 24 cm -3 (r s <0.9), this was found to be 
too soft. Accordingly, a new harder GGA norm conserv- 
ing local pseudopotential was generated using the opti- 
mized pseudopotential method as included in the Opium 
package 4 ^,. This new potential was then used in the 
CPMD code. It has a cutoff radius of r c = 0.5 &b (q c = 15, 
6 Bessel functions) and requires a plane wave cutoff of 
100 Ha to yield convergent results^. 

MD runs were usually performed using only T-point 
sampling of the Brillouin zone. Corrections due to k- 
point sampling can become important in the metallic 
region of hydrogen. Here, contributions of at most 2% 
were determined at randomly chosen snapshots where the 
pressure was reevaluated using 2x2x2 and 4x4x4 
Monkhorst-Pack grids of k-points in abinil— . Effects due 
to finite electron temperatures (Fermi smearing) have not 
been found for the high densities needed to achieve over- 
lapping with the analytical theory. 

As a first comparison between properties of our two- 
fluid model and ab initio DFT-MD simulations, Fig. [5] 
shows ion-ion pair distributions. The results using VASP 



and CPMD deviate due to differences in the number 
of particles considered and the type of pseudopotential 
used. In the present study, all VASP runs have been per- 
formed with 128 electrons and 432 electrons were used in 
CPMD simulations. In addition, the CPMD runs use a 
much harder electron-ion pseudopotential having a core 
radius of r c = 0.5 a^ whereas r c — 0.8 a^ is used in the 
VASP runs. For the densities with r s = 0.9 and r s = 0.7, 
these differences are insignificant. However, both the 
harder pseudopotential and the increase in the number 
of particles / box size are important for the highest den- 
sity with r s =0.5. Here, VASP predicts a first correlation 
peak that is too high which can be traced back to the too 
large core radius of the pseudopotential. Furthermore, a 
box with 128 particles is too small to allow for the com- 
putation of the structure at larger distances as it is seen 
in the tail of the pair distribution. 

The pair distribution functions obtained from solutions 
of the hypernetted chain equations (HNC) are consistent 
with the approximations made to derive the two-fluid 
model. The HNC solver uses a linearly screened ion- 
ion potential and keeps the nonlinear correlations within 
the ion subsystem only. The comparison with the DFT- 
MD results shows that this is an appropriate model for 
the high densities shown (note that better agreement is 
obtained with the CPMD data for r s =0.5). Here, the 
electron-ion coupling strength is sufficiently low to justify 
the linear response formalism. For the lowest density 
with r s = 0.9, there is some deviation in the first slope 
and in the strength of correlation that are a result of the 
linear screening approximation used. The comparison 
with the DFT-MD data demonstrate that the model of 
linear screening is beyond its limit of applicability for this 
density. 



IV. RESULTS FOR THE EOS 

We have already established that our two-fluid model 
merges with the exact T = limiting law at the high 
density site of the EOS. In the following, we will show 
that this model and improved DFT-MD simulations give 
results in agreement with each other over a range of den- 
sities and temperatures. As DFT-MD and further Quan- 
tum Monte Carlo simulations already meet the exact low 
density limit (Dcbyc-Huckel or density/fugacity expan- 
sions), the entire density range can now be described by 
theories in the pysical picture. 

Figurc|3]and Table Q] include a comparison of EOS data 
obtained by the different simulation techniques. CEIMC 
results were taken from recent worked. Unlike earlier 
CEIMC results^, they give almost identical values for 
the pressure as obtained from DFT-MD. The data points 
span almost an order of magnitude in density and cover 
the region of metallic hydrogen which is most difficult to 
describe. This excellent agreement gives large confidence 
in both first principle simulations. As both data sets have 
a similar slope, we can expect that either both or none 
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FIG. 3: (Color Online) Ratio of pressure to ideal pressure 
in hydrogen at T = 10 4 K versus density. Data shown are 
obtained using VASR 1 ^, VASP 1 ^ up to r s = 0.8, CEIMG^, 
Saumon & Chabrier EOS^^, and SESAMES. The two-fluid 
model uses an effective ion-ion potential in RPA. 

merge with the high-density description of the two-fluid 
model. 

Figure [3] also shows a comparison of data obtained by 
ab initio simulations and results of the two-fluid model. 
At the low density side where CEIMC data exist, one 
finds considerable deviations. These differences are due 
to the decoupling of electrons and ions in the two-fluid 
EOS which is not applicable here. 

To test the merging of the simulation data into the 
two-fluid description, we have to consider higher densi- 
ties where the electron-ion coupling is weaker. For this 
task, we turn to DFT-MD simulations. In the case of 
VASP, this is not straightforward as one is not free in 
the choice of the electron-ion pseudopotential. Particu- 
larly, for densities above r s =0.7, the particles are closer 
together than the core radius of the pseudopotential. As 
we have already shown for the ion-ion pair distribution 
(see Fig. [2]), this causes VASP to yield unreliable results 
for r s < 0.7. As a result, the pressure obtained from 
VASP first approaches the results of the two-fluid model 
and then starts to deviate again for high densities. 

The behavior of the VASP data demonstrates the need 
for a harder pseudopotential with a smaller core radius. 
Moreover, finite size effects become important for higher 
densities with larger correlation lengths. Both issues have 
been resolved in CPMD runs using a new pseudopotential 
with a core radius of r c = 0.5 &b and N = 432 electrons 
and protons. The result show a smooth merging with the 
two-fluid model at high densities. 

Figure [3] shows also results from EOS models often 
used in planetary modelling namely the SESAME ta- 
bles and the EOS model constructed by Saumon & 
Chabrier—. The latter perfectly merges into our two- 
fluid model, but deviates from the quantum simulations 



TABLE I: Comparison of the hydrogen EOS as obtained 
by VASP 1 ^^, CEIMG24, CPMD (this work, using an im- 
proved electron-proton pseudopotential) and the two-fluid 
model (last two columns) using two effective proton-proton 
potentials as indicated. 
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VASP 


CPMD 
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0.589 
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at lower densities. The SESAME data, on the other 
hand, disagree with both ab initio simulations and the 
two-fluid model over the whole density range considered. 

V. EXTENSION AND DISCUSSION 

A. Nonlinear Electron-Ion Interactions 

A limitation of the two-fluid model arises from the use 
of first order perturbation theory to describe the response 
of the electrons to the fields created by the ions. One can 
however argue that the good agreement of the two-fluid 
model and the fully nonlinear simulations is a result of 
the fact that quadratic response terms cancel to a large 
extent^. 

An improved treatment including the fully nonlin- 
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FIG. 4: (Color Online) Pressure of dense hydrogen at T = 
3500 K normalized by the ideal pressure. The symbols mark 
simulation data similar to Fig. (CEIMC data are interpo- 
lated for this temperature). The two curves show results of 
the two-fluid model with different ion-ion potentials: the lin- 
early screened Yukawa potential (solid, red line) and the non- 
linear Hulthen potential (red, dash-dotted line). 



ear response (see, e.g., Refsi^— ) is beyond the scope 
of this paper. Here, we estimate nonlinear effects in 
the electron-ion interaction by applying the nonlinear 
Hulthen potentials! 



2 

6 K e 



i 



(13) 



as an ad-hoc model for the effective ion-ion interaction. 

Interestingly, the results of the two-fluid model agree 
rather well with the data from the quantum simulations 
for 3500K in FigfJ] if this nonlinear potential is used. 
In the density range with 1.1 < r s < 0.7, the agree- 
ment is much better than for the case of the linearly 
screened Yukawa potential. Indeed, Potekhin & Chabrier 
include local field corrections in the screening for the ef- 
fective electron-proton interactional. For the parame- 
ters of FiglH our curve using the Hulthen potential and 
the curve according to their model are indistinguishable. 
However, there is no effect due to nonlinear electron-ion 
interactions for temperatures above 5000K and densities 
larger than r s = 0.9. 



B. Quantum Effects in the Ion Subsystem 

For the highest densities considered, quantum effects 
may also become important for the ions. To estimate 
quantum effects on the protons, we first consider a po- 
tential that accounts for quantum diffraction effects of 



order e 2 — 



«P(r) 



Z 2 e 2 ^ 



exp —K e r 



(14) 



r_ _ AjKe 

Xi 2 



exp K e r 



2$ 



1 - a! 



- 1 



Here, n e is the inverse screening length of the electrons, 
Xi = y/ a/ ' (2mukBT) is the thermal wavelength of the ions 
with the reduced ion mass mu = m,i/2, and <&(x) denotes 
the error function^. For large distances, the quantum 
potential (jT4"]) approaches the screened effective ion-ion 
potential (|S|). At the origin, it has the finite value 



lim v 

r->0 



KK 



2^2 



Z 2 e 



X, 



■ exp 



X 2 n 2 



1 - $ 



X,K, 



(15) 

which reflects quantum diffraction effects. 

Moreover, quantum exchange effects can be important 
for high densities. These can be included by adding the 
following exchange potential^ 



£ (r) = iexp(-r 2 /A 2 ) 



(16) 



k B T- 



Z 2 e 2 7r 
4r 



da , 
— $ 

a 



The first term accounts for ideal exchange in an averaged 
way whereas the second term gives the e 2 contribution. 

The quantum potentials fT3|) and (Tpp)) are derived from 
Slater sums, and are exact in the sense of perturba- 
tion theory. They give the correct quantum thermo- 
dynamic functions in the weakly coupled and weakly 
degenerate limit, i.e., for systems with Tu -C 1 and 
nAf = n(2irh/m,ikBT) 3 / 2 <C 1. We use these potentials 
here to estimate for which densities quantum effects start 
to influence the ion properties. It should be emphasized 
that the quantum potentials (|14|) and (|16p are used only 
in the classical method employed to determine the pair 
distribution function gu (here, in the HNC solver). The 
thermodynamic expressions @ and ([¥]) are valid in quan- 
tum as well as in the classical case. Hence, the effective 
ion-ion potential ([5]) must be used without any quantum 
corrections in these formulas. 

Figure [5] shows ion pair distributions obtained using 
the usual screened potential and the quantum pseudo- 
potential. Deviations become obvious for plasmas with 
nAf > 0.6. It can be observed that the quantum effects 
weaken the short range order in the proton system. It is 
however remarkable that these effects only set in at these 
high values of degeneracy. The reason for this behavior is 
found in the strong Coulomb forces at high densities that 
create a large correlation hole. Consequently, the protons 
are too far apart to experience short range diffraction and 
quantum exchange effects. 
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FIG. 5: (Color Online) Ion-ion pair distributions for T= 10 4 K 
and three densities calculated by HNC using the screened po- 
tential © (solid lines, labelled Y) and the screened quantum 
potential with exchange (|14l) plus (|16[) (dash-dotted lines, 
KK+ex). The ion degeneracy is nAf = 0.26, 0.52, 1.05 
for the three densities, respectively. The ion-ion coupling 
is Tu = 99(96), 125(119), 157(144), respectively, where the 
quantum coupling strength is given in brackets. The upper 
panel gives the ratio of the pair distributions with and without 
quantum effects. 



The emerging quantum effects in the ion structure 
yield only small changes in the ion contribution to the 
thermodynamic functions. Wc determine a reduction of 
less than 1% in the correlation contribution to the ion 
pressure for plasmas with nAf < 1. At higher densities, 
the deviations from the classical result can be large, but 
their exact calculation requires PIMC methods. 

The only quantum effect outside our control are zero 
point oscillations of the ions. In strongly coupled liquids, 
caging of the ions occurs and, for the time the cage is sta- 
ble, the ions perform oscillations as in a solid^i. PIMC 
calculations in the solid and fluid phases of a Yukawa sys- 
tem did indeed find contributions due to the zero point 
motion, but for very high densities of n = 10 27 cm"—. 
Moreover, the overall change in energy due to quantum 
oscillations is below 1% for the parameters considered. 
Furthermore, DFT-MD simulations, that do not include 
quantum oscillations of the ions, and CEIMC results, 
that include them, agree well in the considered parameter 
range which again indicates that these effects are small. 
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FIG. 6: (Color Online) Pressure over ideal pressure of a hydro- 
gen plasma with T= 10 4 K for a wide density range covering 
the ideal classical plasma state, the atomic gas, the molec- 
ular gas, and the metallic fluid (from left to right). The 
lines show the fugacity expansion from Ref.— , the two-fluid 
model (this work) using an effective ion-ion potential in RPA, 
and the T = limit. The symbols mark simulation data 
from RPIMO^, CEIMO^, and DFT-MD from this work and 
Refs^. 



C. Hydrogen EOS for Arbitrary Densities 

The two-fluid model, based on the quantum statistical 
Green's function approach for the electron contributions 
and classical HNC methods for the ion properties, is well 
applicable in the density range that is neither covered by 
present first principle simulations nor by the T = limit. 
After examining the high-density part of the EOS in de- 
tail, it is worth noting that theories and simulations in 
the physical picture can be used to cover the entire pa- 
rameter range and show agreement in overlapping regions 
of their applicability. 

Figure [S] demonstrates this fact for an isotherm that 
covers the low density ideal plasma, the atomic fluid, the 
molecular fluid, and the metallic fluid to the border of 
the fully degenerate electron-proton system. Although 
each phase requires well-suited theories or simulations, 
one finds a smooth EOS that is often based on several 
methods and that does not involve any interpolation. 
Such a combined EOS can serve as an excellent basis for 
the modelling of gas planets and other compact objects 
dominated by hydrogen. 



VI. CONCLUSIONS 

We have applied a two-fluid model to investigate the 
hydrogen EOS at high densities. Our approach combines 
the advantages of a quantum theory based on thermody- 
namic Green's functions for the electrons with a classical 
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description of the structure in the ion fluid. Thus, it is 
able to account for electron degeneracy, finite tempera- 
tures, and strong forces between the ions. Its only limita- 
tion is the requirement of weak electron-ion interactions 
which naturally excludes bound states. 

The two-fluid model agrees very well with a number of 
exact limiting laws: i) at low densities it almost exactly 
merges with the Debye-Hiickel law; ii) at high energies 
it practically coincides with the T = limit as long as 
quantum effects on the proton subsystem are negligible; 
and iii) for high plasma temperatures it agrees with a 
perturbation expansion in terms of interaction strength 
for the entire density range. 

In the high density region, we find excellent agreement 
of the two fluid model with first principle simulations 
(DFT-MD). Our DFT-MD simulations agree also well 
with recent CEIMC results. At lower densities (r s >0.7), 
the request of weak electron-ion coupling is not fulfilled 
and one finds increasing deviations between the results 
of the two-fluid model and the quantum simulations. For 
densities with r s < 0.7, the two-fluid model is however 
a reliable and computationally cheap alternative to full 
scale quantum simulations. From the agreement in pres- 
sure and ion structure between simulations and two-fluid 
model, we can conclude that the exchange correlation 
functional used in DFT-MD and our Green's function 
approach give the same electronic contributions for the 



relevant parameters. The two-fluid model also bridges 
the parameter space where neither DFT-MD simulations 
nor the T = law can be applied. This success means that 
hydrogen can confidently be described in large parts of 
the phase diagram by techniques using the physical pic- 
ture. For temperatures of a few electron volts relevant for 
astrophysics and inertial fusion, this includes the entire 
density range. 

Although the two-fluid model was applied here only for 
hydrogen, it is also applicable for any fully ionized system 
with higher ion charge states or very stable inner shell 
configurations with small modifications. The essential 
test is if the electron-ion interaction can be considered 
to be small. Accordingly, systems like the fluid region 
in white dwarfs can be confidently described by the two- 
fluid model presented here as well. 
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